This invention relates to generally to control systems and more particularly to model predictive control systems.
Model predictive control logic has many advantages in controlling practical multivariable systems. A multivariable system has multiple goals (also known as commands) and has multiple effectors that change system dynamics, for achieving those goals. Multivariable systems may also include significant cross-coupling, in which key effectors each significantly effect multiple goals. Therefore, a controller for this system should also be multivariable to decouple the responses. In a decoupled system, when a goal changes, the control algorithm responds by setting multiple effectors so that this one goal changes without significantly changing the other goals. Furthermore, the system cross-coupling is often dynamic, changing in time, where some coupling is persistent, some transitory, some fast, some slow, or it may first couple in one direction, then reverse and couple in the other. Therefore the control algorithm decoupling should also be dynamic, evolving in time to counter the changing nature of the coupling. In addition to de-coupling, the multivariable control logic must also provide the tradition feedback properties of reducing the effects of changes in the system dynamics and the effects of disturbances, all while ensuring the system remains stable, safe, and durable.
There is much control theory, some reduced to practice, supporting the design of feedback control logic for coupled multivariable systems. Most of this theory presumes the system to be controlled has no inequality constraints, referred to as limits here, on any of the system variables or effectors. This is a very significant mismatch between theory and practice. Real and practical systems have substantial and significant limits. Some limits are physical, such as effector ranges, and some are operational, to ensure stability, safety, and durability. In a real sense, if no limits are hit, the system has been over designed. This mismatch has led to great difficulties and ad hoc methods in control law design, including splitting logic into many modes, where there is a separate control mode to cover each permutation of limit or fault condition, and more ad hoc logic to dynamically control the switching between separately designed dynamic modes. Coordinating these modes and ensuring there are no adverse dynamic interactions between modes (e.g., repeatedly jumping back and forth between two modes, sudden jumps, etc.) requires a significant amount of additional ad hoc logic. Ad hoc design makes the control software difficult to design, difficult to maintain and greatly reduces the degree to which the software can be reused on a similar system.
Model predictive control (MPC) is a multivariable control theory that explicitly includes limits and, therefore, provides a good match with practical systems. MPC can also be configured to respond in realtime to changes in the system, such as actuator faults. Thus, MPC provides a formal method of control algorithm design for multivariable systems that can decouple responses, as well as physically possible, even as limits are hit and faults occur. In MPC design, there is no need for ad hoc designs of extra control modes to handle limit conditions and faults. This is expected to provide significant increases in controller software design productivity, make it easier to upgrade and otherwise maintain control software, and make it possible for the control system to accomplish more difficult tasks. The latter include more autonomous operation, more completely integrating system responses, and operating with reduced physical stresses.
The most significant roadblock to implementing model predictive control is the amount of computation required. That is why the initial practice of MPC, which began 20 to 30 years ago in the chemical industry, was for systems with relatively slow dynamics. A chemical plant controller updates actuator settings on the order of once every 5 minutes. These early applications stressed the computer power available at the time. In a gas turbine engine, the update rate must be around 40 times per second, an approximately 12,000-fold difference, and the control function may be more complicated. Even so, MPC may be practical for engines, and other vehicle applications, because computers are now available that are 500 to 1000 times faster. That still leaves a 10 to 30-fold gap, depending on system complexity, between the computation MPC appears to require and computer power available. This invention is critical because it provides numerical techniques for reducing the computational burden of model predictive control by 20 to 100-fold, closing the gap.
Model Predictive Control
Model predictive control theory has been reduced to practice, most notably in the chemical industry. In support of this, model predictive control theory has been described in numerous technical papers and textbooks. It is summarized here. The main objective of this summary is to show that implementing MPC requires a large amount of onboard computation and that approximately 90% of this computation involves numerous solutions of a large matrix equation in realtime. This section shows that when the optimality equations and variables are ordered in a certain way, the matrix equation has a particular sparse, banded, symmetric structure. These properties imply the matrix includes many zeros, the nonzero entries are gathered near the diagonal, and a matrix “square root” exists. This particular form is developed in this section and termed the Large Sparse Matrix Equation (LSME) equation. The square root method exploits this structure in a novel way to greatly reduce computation.
In model predictive control, actuator commands, u, are determined as the result of solving an optimization problem in realtime. Model predictive control is a type of digital controller. Digital controllers are based on a computer that periodically samples operator commands and measured responses. As a consequence of each sample, the digital controller specifies new actuator settings, which is referred to as a controller update. These settings are most often held until the next update. MPC specifies these settings as the first time point of an actuator trajectory. The trajectory specifies u(n), for controller update times numbered with integer n, where n varies from L to N+L−1. That is, for N time points. The current sample time is enumerated L. This actuator trajectory is determined as the unique trajectory that minimizes the performance index, PI, where
                    PI        =                ⁢                                            1              2                        ⁢                          x                              N                +                L                            T                        ⁢                          M                              N                +                L                            T                        ⁢                          M                              N                +                L                                      ⁢                          x                              N                +                L                                              +                                    1              2                        ⁢                                          ∑                                  n                  =                  L                                                  N                  +                  L                  -                  1                                            ⁢                              [                                                                                                    x                        T                                                                                                            u                        T                                                                                            ]                                                                                    ⁢                                            [                                                                                          C                      T                                                                                                                                  D                      T                                                                                  ]                        ⁢                          Q                              T                /                2                                      ⁢                                                                                Q                                          1                      /                      2                                                        ⁡                                      [                                                                                            C                                                                          D                                                                                      ]                                                  ⁡                                  [                                                                                    x                                                                                                            u                                                                              ]                                            n                                +                      [                                                                                f                    n                    T                                                                                        g                    n                    T                                                                        ]                                                          ⁢                                            [                                                                    x                                                                                        u                                                              ]                        n                    +                                    u              T                        ⁢                          R                              T                /                2                                      ⁢                          R                              1                /                2                                      ⁢                          u              n                                          
and where the summation is over the N future time points. xn is a vector of the states of a system dynamic model and un are the actuator commands at each of the N time points. The matrices C and D are coefficients of linearized output equations, which relate system outputs for which there are goals to the dynamic model states and actuator settings. Q and R are lower triangular weighting matrices, and vectors f and g are driving terms that are functions of a desired trajectory for the system output, desired or nominal actuator values, and weighting matrices Q and R. Lower triangular matrix M is a terminal-weighting matrix and may be present for a variety of control theoretic reasons.
The nominal system and, therefore, the optimization must obey the state dynamics model of the systemxn+1=An·xn+Bn·un+bn
which optimization procedures term equality constraints. The optimization must also obey the physical and operational limits linked to theLimited variables yc(n)=Ccn·xn+Dcn·un+an≦Yn
for the N time points. The first equality represents a model of how the limited variables are linked to the state dynamics model. The explicit inclusion of the limit equation in the theory and resulting logic design is a key advantage of MPC.
A MPC designer specifies the control algorithm by determining the state dynamics model, the limit model, and the performance index. Performance index design involves selecting the variables to include in output vector, y, and choosing the weighting matrices Q, R, and M. The control software designer must also provide reusable software for solving the optimization problems in realtime.
The above optimization problem is a quadratic programming problem: a special type of constrained mathematical optimization. Quadratic programming is much more reliable than general purpose constrained optimization, and the associated restrictions are mild. The two major methods for solving quadratic programming (QP) problems on a computer are Active Set (AS) methods and Interior Point (IP) methods. The invention described below will significantly reduce the amount computation needed with either.
In MPC operation, the state vector at the current sample time, xL, is given. It is a function of measurements of the system response up to the controller update time L. Thus, the MPC controller is a feedback controller. The net output of the above optimization is the vector uL, the first time point in the u trajectory, which is the actuator command for the Lth update of the MPC controller.
The same MPC optimization must be performed, with the variables shifted in time, to determine the next actuator command, uL+1. The MPC optimization must be re-solved every controller update period. The next three sections summarize how the above type of optimization is most commonly solved in an analytical fashion.
Adjoint Methods for Absorbing Constraints Into J
In the usual method of solution for the optimal control, the equality constraints (the state dynamics model) are adjoined to the performance index with Lagrange multipliers, pn, and the inequality constraints (the model of system variables with respect to their physical limits) are adjoined with multipliers, mn, to form the augmented performance index as follows.
                    Ja        =                ⁢                                            1              2                        ⁢                          x                              N                +                L                            T                        ⁢                          P                              N                +                L                                      ⁢                          x                              N                +                L                                              +                                    ∑                              n                =                L                                            N                +                L                -                1                                      ⁢                                                                                1                    2                                    ⁡                                      [                                                                                                                        x                            T                                                                                                                                u                            T                                                                                                                ]                                                  ⁢                                [                                                                                                    C                        T                                                                                                                                                D                        T                                                                                            ]                            ⁢                              Q                                  T                  /                  2                                            ⁢                              Q                                  1                  /                  2                                                                                                    ⁢                                                            [                                                                            C                                                              D                                                                      ]                            ⁡                              [                                                                            x                                                                                                  u                                                                      ]                                      n                    +                                                    [                                                                                                    f                        n                        T                                                                                                            g                        n                        T                                                                                            ]                            ⁡                              [                                                                            x                                                                                                  u                                                                      ]                                      n                    +                                    u              T                        ⁢                          R                              T                /                2                                      ⁢                          R                              1                /                2                                      ⁢            u                    +                                                ⁢                                            p                              n                +                1                            T                        ⁡                          (                                                                    A                    n                                    ·                                      x                    n                                                  +                                                      B                    n                                    ·                                      u                    n                                                  +                                  b                  n                                -                                  x                                      n                    +                    1                                                              )                                +                                                ⁢                              m            n            T                    ⁡                      (                                                            Cc                  n                                ·                                  x                  n                                            +                                                Dc                  n                                ·                                  u                  n                                            +                              a                n                            -                              Y                n                                      )                              and the terminal weight PN=MNTMN.
In the Lagrange method, the adjoined constraint terms are formulated so that the constraints are satisfied when the augmented performance index is optimized as if there were no constraints. The adjoined terms also have zero value at the optimal. Thus adjoining the constraints does not change the value of optimal J. For inequality constraints the adjoint method also requires that for each limit, at each time n, either the corresponding m=0 or the inequality constraint is on its limit.
Karush-Kuhn-Tucker (KKT) Conditions
A set of algebraic equations whose solution specify the optimal system trajectory are known as the (KKT) conditions. For the above MPC optimization, it is known that the solution includes the following KKT conditions for each controller update between time L and time N+L−1:
                                                                        [                                                                                                    C                        T                                                                                                                                                D                        T                                                                                            ]                            ⁢                              Q                                  T                  /                  2                                            ⁢                                                                                          Q                                              1                        /                        2                                                              ⁡                                          [                                                                                                    C                                                                                D                                                                                              ]                                                        ⁡                                      [                                                                                            x                                                                                                                      u                                                                                      ]                                                  n                                      +                                          R                                  T                  /                  2                                            ⁢                                                R                                      1                    /                    2                                                  ·                                  u                  n                                                      +                                                                          [                                                                    f                                                                                        g                                                              ]                        +                                          [                                                                                                    A                        T                                                                                                                                                B                        T                                                                                            ]                            ⁢                              p                n                                      -                                          [                                                                            I                                                                                                  0                                                                      ]                            ⁢                              p                                  n                  +                  1                                                      +                                          [                                                                                                    Cc                        t                                                                                                                                                Dc                        t                                                                                            ]                            ⁢                              m                n                                                          =    0                                A          n                ·                  x          n                    +                        B          n                ·                  u          n                    +              b        n            -              x                  n          +          1                      =    0                                Cc          n                ·                  x          n                    +                        Dc          n                ·                  u          n                    +              a        n            +              t        n            -              Y        n              =    0                      t        n            ·              m        n              =    0  the split boundary conditions:xL=xfeedback and pN+L=MN+LTMN+LxN+L 
and the additional conditions that all t's and m's must be greater than zero. If these equations are solved, the actuator commands, u(n), that are part of the solution, used as inputs to state equation model of the system dynamics, will minimize the performance index, while meeting the limits.
General Large Sparse Matrix Equation for Quadratic Programming
The KKT conditions for the N time points can be collected into the following large sparse matrix equation (LSME). This must be repeatedly solved when solving the MPC quadratic programming problem, using either Interior Point (IP) or Active Set (AS) methods. The LSME has the form:
            [                                                                  -                                  T                  T                                            ⁢              T                                            J                                                              J              T                                            H                              ]        ·          [                                    m                                                z                              ]        =      [                            K                                      f                      ]  
where z is a vector containing the state, x, control, u, and state equation adjoint variable, p, for each time point, group by time. m is the adjoint variables for the inequality constraints, grouped by time. f and K are vectors. H and J are banded matrices, where
      H    ≡          [                                                                                    D                  0                  T                                ⁢                                  Q                  0                                ⁢                                  D                  0                                            +                              R                0                                                                        B              0              T                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          B              0                                            0                                              -              I                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0                                              -              I                                                                          C                1                T                            ⁢                              Q                1                            ⁢                              C                1                                                                                        C                1                T                            ⁢                              Q                1                            ⁢                              D                1                                                                        A              1              T                                                                                                                                                                                                                                                                                                                                                                                                                        0                                                              D                1                T                            ⁢                              Q                1                            ⁢                              C                1                                                                                                          D                  1                  T                                ⁢                                  Q                  1                                ⁢                                  D                  1                                            +                              R                1                                                                        B              1              T                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                A              1                                                          B              1                                            0                                              -              I                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0                                              -              I                                                                          C                2                T                            ⁢                              Q                2                            ⁢                              C                2                                                                                        C                2                T                            ⁢                              Q                2                            ⁢                              D                2                                                                        A              2              T                                                                                                                                                                                                                                                                                                                                                                                                                        ⋰                                                                                                        ⋰              ⁢                                                                                                                                                  ⋰                                                                                                                                                                                                                                                                                                                                                                                                            0                                                              D                                  N                  -                  1                                T                            ⁢                              Q                                  N                  -                  1                                            ⁢                              C                                  N                  -                  1                                                                                                                          D                                      N                    -                    1                                    T                                ⁢                                  Q                                      N                    -                    1                                                  ⁢                                  D                                      N                    -                    1                                                              +                              R                                  N                  -                  1                                                                                        B                              N                -                1                            T                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                A                              N                -                1                                                                        B                              N                -                1                                                          0                                              -              I                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0                                              -              I                                                          Q              N                                          ]                  J      =              [                                                            Dc                1                                                    0                                                      Cc                1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            Dc                2                                                    0                                                      Cc                2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    ⋰                                      ⋰                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              Dc                N                                                    0                                                      Cc                N                                                    ]              ,          z      =              [                                                            u                0                                                                                        p                1                                                                                        x                1                                                                                        u                1                                                                                        p                2                                                                                        x                2                                                                                        u                2                                                                        ⋮                                                                          x                                  N                  -                  1                                                                                                        u                                  N                  -                  1                                                                                                        p                N                                                                                        x                N                                                    ]              ,                  and        ⁢                                  ⁢        m            =              [                                                                                                                                                                                                      m                            1                                                                                                                                                                            m                            2                                                                                                                                                                                    ⋮                                                                                                                          m                N                                                    ]            
T is a diagonal matrix and, the matrices are shown for the controller update at time L=0. Matrix H is also symmetric, which implies there is a square root matrix, Hr, such that Hr′*Hr=H. The LSME has the same form at other time points but with different parameter values. The LSME is common to both of the major quadratic program solvers, except that the value of T and the right hand side (RHS) vector depend on the quadratic program solver used. With either method the LSME must be solved repeatedly (10 to several hundred times, depending on system complexity) for each controller update.
This form of the LSME is novel and special. It is achieved by ordering the optimality equations (rows of the LSME) and the variables in z in a special way. Then the resulting LSME matrix is sparse, banded, and symmetric.
Interior Point Methods
In the Interior Point method for solving quadratic programs, T is diagonal. Each element is the square root of the product of the corresponding slack variable, t, and adjoint variable, m, for each time point. J is sparse banded matrix, and includes all the limit equations. For the most popular family of interior point methods, T and the right hand side vector vary twice each iteration. The right hand side vector is determined from the solution of the previous iterate. H and J remain constant for the iterations at the same controller update period. The solution for the optimal actuator settings will require 10 to several hundred iterations, depending on problem complexity. For more explanation, see Ref. 1, especially their case where the slack variable has been eliminated.
Active Set Methods
The Active Set method searches for which limits are active. Only some of the inequality constraints, Jc*z<=Ymax, will be active. Each iteration of an Active Set qpsolver “guesses” which of the constraints are active. Let the vector “active” contain the indices of the active constraints. When active, the inequalities become equalities, and Jc(active,:)*z=Ymax(active), and the system is running on the active limits. Compared to the LSME, J=Jc(active,:), K=Ymax(active), and T=0. In a given controller update period, a number of iterations of “active” will be evaluated in a search for the correct active set. The search strategy determines a new iteration of “active” based on the results of the previous iterate's solution of the LSME. H will remain constant for the iterations for a given controller update.
The key observation here is that some form of the LSME must be solved numerous times to solve the MPC quadratic program each controller update, whether interior point or active set solver is used. For practical controllers, the LSME is solved 10 to 100 times per controller update, depending on N and other factors.